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1 Introduction 

Experiments aimed at the first direct detection of gravitational waves are now 
underway with a network of laser-interferometric gravitational-wave (GW) 
detectors [l][2l[3] . At the completion of the current science run the LIGO and 
Virgo detectors will be upgraded to second generation detectors, with an order 
of magnitude improvement in sensitivity, and Advanced LIGO is due to go 
online in 2014. It is hoped that the second-generation detectors will yield 
enough observations to begin in earnest the field of GW astronomy; see [1] for 
an overview of the astrophysics and cosmology potential of GW observations. 
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With the goal of achieving a further factor of ten improvement in sensi- 
tivity, combined with an extension of the detector bandwidth to the range 
IHz to lOkHz, a design study for a third generation detector has been sup- 
ported within the European FP7 framework. If these design goals are met, the 
Einstein Telescope (ET) would be sensitive to a volume of the universe one 
million times larger than current ground-based detectors, and allow a precision 
probe of gravitational-wave sources throughout the universe to cosmological 
distances. 

Extracting the most science from ET observations will require accurate 
theoretical knowledge of GW sources — far more accurate than that required 
for first- and second-generation detectors. One important aspect of source 
modelling is the numerical solution of Einstein's equations for the inspiral and 
merger of compact bodies (black holes and neutron stars) and the numerical 
simulation of the physics of single compact bodies (neutron star oscillations, 
instabilities and collapse, and core-collapse supernovae amongst others). 

NR simulations have advanced significantly in the last few years, partic- 
ularly in describing the last orbits and merger of systems of two black holes 
(see OISllT] for the original breakthrough results, and [8] for an overview of 
the status of current simulations focused on GW detection). The mode-ling 
of neutron-star binaries has also made significant progress from the first sim- 
ulations of the Tokyo group ([9]) to the current status as recently reviewed 

by m- 

In this paper we will discuss some of the issues that numerical-relativity 
simulations face in producing the results needed for ET data analysis and 
astrophysics applications. 

For black-hole binary simulations, the situation is easy to state: since the 
total mass of the system provides an overall scale to the numerical results, 
the increased bandwidth of ET does not change the physical nature of the 
simulations that need to be produced; the black-hole binary parameter space 
remains the same. However, it docs change the accuracy requirement of the 
simulations, and of the model of the GW parameter-space that can be pro- 
duced from them. 

For matter simulations, the situation is quite different. The scales are set 
by the physics included in the model. The size of the parameter space is 
dominated by the unknown physics which it is hoped that GW observations 
will constrain. In order to make simulations practical, only the "essential" 
physics will be included in the model. The determination of "essential" is 
balanced between current knowledge and expectations of the key physics, and 
what can practically be simulated and detected. Thus the model, resulting 
parameter space and accuracy requirements will change with the increased 
bandwidth of ET. 

In Section [2] the status of vacuum simulations of binary black holes is 
discussed, including a discussion of the waveform accuracy of current and near- 
future simulations, the sampling of the parameter space, the construction of 
template banks and possible future improvements in numerical techniques. 
Section [3] discusses the situation with matter simulations, assuming that the 
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physical model is "added on" to the spacetime simulation and GWs extracted 
using the same techniques. Section [4] summarises the likely status of NR over 
the next decade. 



2 Black-hole binaries 

Fully general-relativistic simulations of the dynamics of black- hole binary sys- 
tems during their last orbits and merger have been possible since 2005 [SI 
[6ll7]. In these simulations initial data are prescribed for a slice of spacetime 
that contains two black holes, usually in orbit, and in many cases with initial 
parameters chosen, or tuned, such that the black holes are following non- 
eccentric quasi-circular orbits [TT | [T2 P T3]. The initial data must also satisfy 
four constraint equations in order to be part of a valid solution of Einstein's 
equations. The data are advanced to future time slices via a set of evolution 
equations; the specific equations depend on the choice of decomposition of the 
four-dimensional Einstein equations to 3-1-1 (space-ftime) dimensions, and the 
coordinates and time-slicing at successive evolution steps depend on the choice 
of gauge conditions. We are in principle free to choose the gauge conditions as 
we wish (as a result of the coordinate invariance of Einstein's equations), but 
in practice are limited to choices that lead to sufhciently stable and accurate 
numerical simulations. The issues involved in finding a numerically well-posed 
and stable formulation of Einstein's equations, a convenient geometrical rep- 
resentation of a black-hole spacetime, construction of black-hole-binary initial 
data, and numerically and geometrically well-behaved gauge conditions, are 
discussed in the textbook [Mj , and in the review articles [IS l fTH l ITT] . 

The parameter space of inspiraling black-hole binaries consists of ten pa- 
rameters: the total mass and mass ratio of the system (or alternatively the 
individual masses of the two black holes), the spin vector of each black hole, 
the eccentricity of the system, and the initial separation of the binary (or al- 
ternatively the initial phase of the GW signal). Since the total mass provides 
an overall scale for the solution, this can be removed from the parameter space 
of necessary numerical waveforms. If the numerical late-inspiral-plus-merger 
waveform is combined with an analytic inspiral waveform (for example, cal- 
culated from the post-Newtonian approximation) , the initial separation of the 
binary can be made arbitrarily large, and yet one more parameter can be re- 
moved from the parameter space. Assuming that the process of constructing 
a "hybrid" waveform is sufficiently robust (and this remains a topic of current 
research), we are left with a total of eight parameters. 

To date late-inspiral-merger-ringdown simulations of black-hole binaries 
have been performed for systems with mass ratios up to 1:10 (although most 
do not extend beyond 1:4), and systems with a variety of spin magnitudes 
and orientations, and for several choices of initial eccentricity. A periodically 
updated summary of simulations that include at least ten GW cycles before 
merger is given in the online version of the review [8]. Note that although a 
large number of simulations have been performed, these only account for an 
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extremely sparse sampling of the full parameter space, and do not yet include 
any simulations of mass ratios beyond q = mxjmi = 10, or of spins above 
a/m ~ 0.92. Efforts to use numerical-relativity simulations to make study the 
astrophysics of electromagnetic counterparts have also begun [TBlfT^ . 

The issues for black-hole-binary simulations, as they relate to GW astron- 
omy, are: (1) determining and achieving sufficient accuracy to extract the max- 
imum physical information from GW observations, which includes producing 
simulations that include enough cycles that subsequent analytic-numerical hy- 
brid waveforms meet these accuracy requirements, (2) producing simulations 
of a sufficiently dense sampling of the parameter space, (3) developing methods 
and codes that are efficient enough to achieve these goals. 

In distinguishing the requirements on numerical waveforms between ET 
and other detectors, for example Advanced LIGO, the key difference is that 
since ET is expected to be roughly an order of magnitude more sensitive than 
Advanced LIGO, the waveform accuracy will also need to be higher to allow 
the most science to be extracted from observations. The question for numerical 
relativity then is what those accuracy requirements are, and whether they can 
be achieved before ET is completed. 

2.1 Waveform accuracy 

In the last few years questions of waveform accuracy requirements have begun 
to be addressed [20], and it seems that current methods and simulations are 
sufficiently accurate for the current generation of ground-based detectors, al- 
though a systematic study has been performed only for the case of equal-mass 
nonspinning binaries |21| . In that work it was shown that, within certain rea- 
sonable caveats, it will not be possible to use GW observations to distinguish 
between the most accurate numerical waveforms used in the study unless the 
signal-to-noise ratio (SNR) is above about 25. 

SNRs of that magnitude are expected to be rare from the Initial and En- 
hanced LIGO and Virgo detectors. However, the Einstein Telescope is expected 
to be an order of magnitude more sensitive, and SNRs in the range 100-1000 
may be typical. We first consider whether numerical simulations will be able 
to match the accuracy requirements of ET by the time of its construction. 

Waveform accuracy requirements are discussed in some detail in [20| . One 
simple criterion for waveforms to be sufficiently accurate for all parameter 
estimation purposes is that the average error in the amplitude and phase 
(appropriately weighted by the power spectrum of the detector noise [20] ) 
satisfy 

S^^'^+SxK^, (1) 

where p is the SNR, and 5x is the fractional error in the amplitude and 54) 
is the error in the phase. If we assume that most black-hole-binary signals of 
interest for ET will have an SNR of less than 1000, then we require that the 
average phase and fractional amplitude errors be of the order of 10"'^. For the 
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waveforms described in [22] (as an example of a state-of-the-art simulation), 
the maximum phase and amplitude errors are roughly an order of magnitude 
larger than the average errors; we expect that the detector-noise-weighted 
average errors will behave similarly. 

As such, the most accurate numerical simulations to date [23] have total 
phase and amplitude errors of 0.02 radian and 0.3% respectively (see the table 
in [21] for a convenient comparison of errors between codes), and so we expect 
that the average errors are already within the accuracy requirements of ET. 
The simulation that deals most convincingly with the wave extraction |22j 
gives average phase and amplitude errors also within the approximate 10"'^ 
requirement provided here. 

So it appears that state-of-the-art numerical simulations are already at, 
or close to, the accuracy requirements for ET science. Assuming that Moore's 
Law holds, then the increase in computing power by the time ET is operational 
will be 2^"/^-^ w 100. Current codes typically employ spatial derivatives that 
converge to fourth-order or better with respect to the spatial resolution, and 
if we assume fourth-order convergence in future numerical results, we will see 
an improvement of two orders of magnitude over current results in the next 
decade. This is a conservative estimate (it assumes that only the computers, 
and not the codes, will improve!), and serves only to demonstrate that accuracy 
of numerical methods will not be the bottleneck in providing waveforms for 
ET science. 

(It is true that this simple analysis does not take into account higher modes, 
for which comparable accuracy is harder to obtain, but for which high accuracy 
will be needed for parameter estimation [24l[25ll26l[27] . However, even the sub- 
dominant modes are well-resolved in the most accurate recent simulations |281 
[29]. 

As we will see in the next sections, the greater challenge is in sampling the 
black-hole-binary parameter space. 

2.2 Phenomenological modelling of black-hole-binary waveforms 

Individual black-hole-binary simulations can only provide a discrete sampling 
of the parameter space. The practical construction of template banks for 
searches, and the use of theoretical waveforms for parameter estimation and 
other follow-up studies, requires waveforms for any values of the black-hole 
parameters. As such, the natural approach is to construct an analytic model 
of the waveforms, and to use the numerically generated waveforms as input 
to calibrate the free parameters of the model. Such models also aim to extend 
the numerical waveforms to an arbitrarily large number of GW cycles before 
merger (or, put another way, extend them to arbitrarily low GW frequency). 

To date there have been two approaches to this problem. One is to start 
with an effective-one-body (EOB) model of inspiral-merger-ringdown wave- 
forms (which in some forms pre-date the existence of full numerical wave- 
forms [30 l [3T ] [32 p 33 y 34j ) . and to both calibrate unknown EOB coefficients to 
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numerical results, and to introduce extra "flexibility parameters" that allow 
yet greater fidelity between the EOB model and the full GR results. EOB-based 
models for nonspinning waveforms have been pursued by several groups |351 
[551[571[551l5MiniHTlH^ , and have reached an impressive level of faithfulness to 
the full GR waveforms [3i|4Tl|42]. 

The second approach has been to construct a purely phenomenological 
ansatz motivated by post-Newtonian (PN) predictions for the inspiral, ap- 
proximate models of the ringdown, and the observed behaviour of the signal 
in the intermediate plunge and merger phase, as seen for example in |43j . This 
model, in which each ingredient is motivated by known physics, is extremely 
flexible, and has been extended from original work on nonspinning binaries |441 
145.46 to binaries with nonprcccssing spins [47j and has been shown to also 
be applicable to neutron-star binaries |48| . This procedure is limited only by 
available numerical simulations, although it remains to be seen how easily a 
"simple" ansatz can be developed for general precessing-spin binaries. The 
nonspinning waveform model is constructed from ten mass-ratio-dependent 
parameters, which can in turn be modelled by second-order polynomials, there- 
fore requiring at least three simulations to define the model. The EOB models, 
by contrast, can be calibrated by only one simulation. It remains to be seen 
how well this economy of calibration extends to the spinning-binary cases. 

Both the EOB and pure phenomenological approaches require that the 
waveform model for the GW cycles before the beginning of the NR wave- 
form (i.e., the PN or EOB model) is accurate. This requires comparison of 
NR waveforms with their PN counterparts in the regions where they overlap. 
NR-PN comparisons have been performed for equal-mass nonspinning bina- 
ries [49 l [50 l [5T | [52] , equal-mass binaries with non-precessing spins [53] , eccentric 
binaries [54j , and one configuration of an unequal- mass binary with precessing 
spins [13] . These results suggest that PN methods perform well up to the point 
where NR simulations begin. Detailed studies of the accuracy of NR-PN hy- 
brid waveforms based on these simulations remain to be made, but first studies 
of the accuracy of equal-mass hybrids at least suggest that they will fulfil the 
detection requirements of current ground-based detectors for masses as low as 
10 Mq [SSj. Complementary results for studies of pure PN inspiral waveforms 
suggest that they in turn will be sufficient for masses up to 12 Mq [56]. How- 
ever, the wider sensitivity band of ET may require that the hybrid waveforms 
include more accurate PN ingredients, or longer numerical simulations. This 
is a topic that deserves further study in the coming years. 



2.3 Sampling of the parameter space 

Whichever phenomenological approach turns out to be most practical, it is 
clear that simulations covering a fine sampling of the black-hole-binary pa- 
rameter space will be necessary. An exhaustive sampling of the parameter 
space is usually presented as a monumental computational challenge. Let us 
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consider the scope of this challenge, and the degree to which it may be met 
by the time ET might be commissioned. 

There are eight parameters in the black-hole-binary parameter space. A re- 
cent phenomenological family of inspiral-merger-ringdown waveforms |47| was 
able to faithfully model binaries with non-precessing spins using an analytic 
ansatz in which the coefhcients were expressed as third-order polynomials in 
the physical parameters of the model (in this case the mass ratio and a param- 
eter that can be approximately considered as the total spin of the binary). For 
simplicity, let us take the relatively optimistic view that third-order polyno- 
mials in each of the parameters will be sufhcient to model the full parameter 
space, and that we therefore require at least four data points in each direction 
to robustly constrain our model. This means that we will require 4^ « 65, 000 
simulations. This would indeed be a monumental undertaking! 

However, this does not take into account symmetries and degeneracies in 
the parameter space. For example, the work in |47j successfully parameter- 
izes non-precessing binaries by only the total spin of the binary — with the 
consequence that if the individual spins aren't necessary to model the wave- 
form, then they will also be difficult to distinguish. This reduces the amount 
of information that can be obtained from a detection based on this waveform 
model, but also reduces the computational cost of producing the model. 

As a second example, the works of [571158] exploit the symmetries of black- 
hole-binary systems to argue that certain quantities in black-hole mergers (the 
mass, spin and recoil velocity of the final black hole) can be described by a 
model that requires only 26 simulations. Their model deals with only quasi- 
circular inspiral (so the parameter space is reduced to seven, not eight, param- 
eters), and uses a second-order-polynomial ansatz, but nonetheless requires far 
fewer simulations than a naive estimation of 3^ w 2200! 

These early attempts to characterize black-hole-binary systems suggest 
that a full mapping of the parameter space may be possible within the next 
decade, if not the next few years. Before becoming overly optimistic, how- 
ever, we should bear in mind that the error analyses of many current simu- 
lations deal only with the dominant harmonic (the £ = 2, m = 2 mode) of 
the waveform, and to date it is only this mode that has been included in 
analytic/phenomenological models of black-hole-binary waveforms. Also, two 
important areas of the parameter space have not yet been accessed by current 
simulations. 

One is black holes with high spin. All published long simulations to date 
have used black- hole initial data that are conformally flat, an assumption that 
simplifies the construction of the data, but limits the spin that can be mod- 
elled. The behaviour of black-hole-binary systems with near-extremal spins 
(which may be the most common astrophysically |59 p 60 ( l6T]). have there- 
fore not yet been studied in detail. Proposals have been made to construct 
non-conformally-flat initial data for both excision-based |621I63| and moving- 
puncture [64] codes, although only the excision data is sufficiently developed 
for use in full binary simulations. With workable proposals in place for high- 
spin data, it is likely that high-spin binaries will be simulated in detail in the 
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next few years, and this region of the parameter space wiU be wch modelled 
by the time of ET. 

The situation is different for large mass ratios. In current simulations, the 
numerical resolution requirements of the simulation are defined by the smallest 
black hole, while the overall scale is set by the total mass. This means that 
the higher the mass ratio {q = mi/m,2 > 1), the more effective resolution 
is required, and the more computationally expensive the simulation. Long, 
accurate simulations suitable for use in GW applications have been performed 
up to g = 6 [39], and a short simulation of 2.5 orbits [65] exists for q = 10. 
Since the simulation cost (both in simulation time and memory usage) scales 
at best linearly with the mass ratio, and more realistically quadratically, long 
accurate simulations of mass ratios g ~ 10 are at the very limit of current 
methods and resources, and q ~ 100 are out of the question. An improvement 
by a factor of 100 in computer power by the time of ET suggests that mass 
ratios q ^ 30 may be possible by that time, but q ^ 100 will require new 
techniques; some possibilities are described in Section \2A[ 

It may ultimately not be necessary to simulate systems with g > 10 in full 
general relativity to accurately model the full parameter space. For example, 
the work in [IT] includes simulations up to only <; = 3, but models all mass ra- 
tios by explicitly incorporating an analytically calculated extreme-mass-ratio 
limit. The effective-one-body (EOB) procedure also includes by construction 
the extreme-mass-ratio limit [30]. However, in order to verify that the high 
mass ratios are indeed being faithfully modelled, it will be necessary to per- 
form at least a few simulations at very high mass ratios. Even if only a few 
such simulations are required to verify the fidelity of a given model (or set of 
models), new techniques will be necessary to produce those simulations. 

2.4 Numerical techniques and room for improvement 

For codes that employ the moving-puncture method, numerical methods have 
not changed significantly since the first simulations in 2005, but the accuracy 
and length of the waveforms has improved significantly. The first waveforms 
covered only a few GW cycles before merger, while the most recent simulations 
cover up to 30 cycles before merger. The first long simulations that allowed 
comparison with post-Newtonian (FN) results had an amplitude uncertainty 
of around 10% and an accumulated phase error of around 1 radian before 
merger [49]. The most recent simulations can achieve an amplitude error of 
0.3% and an accumulated phase error through inspiral, merger and ringdown 
of 0.02rad [21]. 

Typically mesh-refinement is used to resolve both the small-scale features 
near the black holes, the medium-scale physics of the GW signal, and the large 
scales necessary to push the outer boundary of the computational domain far 
from the source. However, improvements in general code efficiency and the 
introduction of higher-order spatial finite-differencing (first sixth-order in |661 
150] . and later eighth-order in [67U13j ) have allowed impressive improvements 
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in both code speed, memory efficiency, and numerical accuracy. More recent 
results by two groups [68j[29] have demonstrated the efficacy of multipatch 
methods, which allow for a far more efficient gridding of the wave extraction 
region, which may lead to yet greater accuracy (of the wave extraction) and 
computational efficiency in future simulations. The results of Reisswig et. al. 
have further achieved Cauchy-Characteristic Extraction (CCE), in which the 
waves extracted at some finite spatial distance from the source are then evolved 
to null infinity, which is where GW signals can be unambiguously defined [22l 
These results indeed suggest that this method produces results that are 
gauge invariant. 

Simulations with variants of the generalized-harmonic formulation have 
been performed by Pretorius [70ll5l[7T| and by the Caltech-CorncU group [72l 
[52] . The latter's SpEC code makes use of a multi-domain spectral method, 
which allows a high degree of accuracy and, like the multi-patch methods 
described above, allow high accuracy with minimal computational resources 
in the wave extraction region. 

These codes represent the current state of the art. However, far greater 
accuracy and efficiency will be needed to meet the science goals of ET. Is 
there room for improvement? 

The answer is certainly Yes. A number of options are yet to be fully ex- 
plored. One that has already received some attention is the use of mixed 
implicit-explicit (IMEX) time integration to allow far larger evolution time 
steps. A first exploration of this approach has already been made in the con- 
text of scalar fields on a curved background [73]. Another approach with a 
related potential benefit is the use of symplectic integrators [TiUTSirfSlTfT] . In 
this approach, conserved quantities in the physical system are maintained by 
construction, and this in general allows for far larger time steps. The canonical 
example is the evolution of binaries in Newtonian gravity, where dramatic im- 
provements in accuracy are possible. General relativistic systems are of course 
far more complex, and there is the added difficulty that energy can be dissi- 
pated as gravitational radiation, but these issues do not preclude the use of 
symplectic methods, and some research on this topic is already underway. 

A further improvement in both accuracy and memory efficiency could be 
achieved by the use of hyperboloidal slices of the spacetime. Here the slices are 
spacelike near the source, but asymptotically null, meaning that they reach null 
infinity on a finite numerical domain. Since the GW phase is constant on a null 
slice, only a finite number of GW cycles need be resolved between the source 
and the outer boundary. This approach would therefore save computational re- 
sources in the wave zone, while at the same time allowing clean GW extraction 
at null infinity. Although hyperboloidal initial data for such systems has al- 
ready been proposed [78[|79] . stable simulations are more challenging. Progress 
has however been made in spherical symmetry and axisymmetry [80] , and hy- 
perboloidal simulations of black-hole binaries may be possible in the next few 
years. 

In the case of moving-puncture simulations, improvements may also be 
possible by better exploiting the natural geometry of black-hole slices in this 
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approach. The method was developed for data that describe the black holes as 
topological wormholes, but it was later found that the slices quickly asymptote 
to cylinders that end at the black- hole throat (also called "trumpets" ) |8i p 82 1 
[53]. It may be possible to modify the method to exploit more explicitly this 
trumpet form of the data. In particular, if the original concept of the puncture 
approach can be resurrected — where the singular nature of the data near the 
black holes is described analytically and only small deviations from this form 
(due to the presence of the second black hole) are described numerically — 
then it may be possible to achieve dramatic improvements in code accuracy 
and efficiency. 

These are a number of avenues for improvement that are already being 
explored, and have the potential by themselves to make possible large numbers 
of long, accurate simulations possible — but of course we expect that by 2020 
there will be many other ideas, and it is quite conceivable that by then the 
standard methods of the field will be radically different from those today. 

3 Matter 

In contrast to simulations of binary black holes, many calculations of gravita- 
tional waves have not used full general relativity (for example most calcula- 
tions of core-collapse as reviewed by e.g. [51] and some binary simulations such 
as [55] )• This is unsurprising; for scenarios such as core-collapse supernovae 
the effect of including the key physics in the model is more important than a 
perfect model of gravity. However, with the recent successes in evolving the 
spacetime as detailed in section [2] the use of full general relativity in calcu- 
lations of gravitational waves is likely to become standard on the timescales 
important for ET. 

In addition we note that as the matter terms do not affect the principle part 
of the spacetime evolution equations there is only one potential complication 
for the spacetime evolution by adding matter. In vacuum evolutions, away 
from singularities, all quantities remain smooth. As the matter quantities may 
become discontinuous, a simulation that includes such shocks will have an 
effect on the spacetime. It is currently unclear what the impact is in GW 
simulations, as the only work to look at this ([55]) is in 1 -I- 1 dimensions. 
As no effect has been seen in any current simulation we will assume that no 
fundamental problem will affect matter simulations relevant for ET. 

From these basic points we will assume that all current GW simulations 
including matter will be able to benefit from the theoretical and numerical 
advances for spacetime and gauge evolution and wave extraction reviewed or 
proposed in section[2J Therefore the bottleneck in computing accurate, realistic 
GW signals will be in the physical model of the matter that can practically 
be evolved, and the accuracy of the numerical methods used to evolve it. 

The parameter space for matter simulations is considerably larger than 
for binary black holes. Firstly the underlying physics of the problem sets the 
mass scale which cannot be removed. Secondly there are a wide variety of 
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GW sources; as well as binary systems there are discs and single star collapse. 
But most importantly there is the uncertainty in the physics that the model 
must include to be realistic (e.g., neutrino transport plays a crucial role in core 
collapse but is much less important in a binary inspiral), and in the details of 
that physics (such as the precise form of the equation of state (EOS) needed 
to close the perfect fluid equations). 

The approach so far has been to focus on limited areas in the parameter 
space to determine which elements of the physical model affect the qualitative 
form of the GWs, and if the GW signal is relatively robust to determine what 
qualitative effects minor changes in the parameters make. The key additional 
parameters in inspirals are currently believed to be the bulk EOS and the 
magnetic field; in core collapse neutrinos and composition are more important 
than magnetic fields. Additional effects such as superfluidity, an elastic crust, 
diffusion, dissipation, or other non-ideal effects may play a role. At present 
even the simplest implementation of these effects is ongoing research. 

The construction of initial data for binary systems is similar to that for 
binary black holes and, as reviewed by [TSlllOj. similar issues of the accuracy of 
the conformal flatness approximation and how to produce quasi-equilibrium 
initial data arise. For single star collapse either axi-stationary solutions are 
perturbed or a pre-calculated model is used, with similar concerns about the 
physical accuracy of the initial conditions. In most simulations small parameter 
studies show that the qualitative behaviour of the gravitational waves is robust; 
one fundamental problem will be discussed later. 

The key question for GWs and in particular for ET is whether numerical 
relativity can compute sufficiently accurate signals. As this crucially depends 
on the physics included in the model to be simulated we will discuss the 
different input physics separately. 

3.1 Hydrodynamics 

Due to the scale and temperature of most interesting GW sources the bulk 
of the matter can be reasonably modelled as a fluid, and in highly dynamical 
situations such as a binary merger the bulk ideal fluid effects, modelled by the 
relativistic Euler equations, will dominate. These equations are typically writ- 
ten in conservation law form [87j as they admit discontinuous (weak) solutions 
for which the conservation law form gives a unique entropy satisfying solution. 
The microphysics under consideration is determined by the EOS which closes 
the system of equations and is typically given as a relation between the pres- 
sure and the independent thermodynamic and composition variables. 

The importance of the EOS is obvious, as has been illustrated by a number 
of simulations (as reviewed by [10], but see in particular [85 ] I88 | I89] ) . Large 
and small scale effects both in terms of the stiffness and the effects of heat in 
simplified form have been covered. However current uncertainty in the EOS is 
large, and at present the approaches followed are to restrict to the "best" EOS 
currently available or to use a simple parametrized EOS ([9^) to best fit the 
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likely candidates. The effectiveness of the parametrized EOS combined with 
numerical simulations has been illustrated by [91] , although it should be noted 
that the accuracy of the simulations discussed there will not be sufficient for 
ET. 

Numerical solutions of the relativistic Euler equations use techniques from 
the standard Newtonian Euler equations. In particular it is currently typical 
to use High Resolution Shock Capturing (HRSC) methods (see e.g. [921193] 
and references therein for reviews). These ensure conservation is preserved 
(necessary to ensure the correct weak solution is obtained if a shock is present) 
and ensure that the solution converges near discontinuities. These techniques 
are nonlinear in the sense that the information used at any given point depends 
on the data in the neighbourhood of that point; this is required to bypass the 
conditions of Godunov's theorem ([94]) and construct numerical methods that 
give correct solutions that converge faster than first order. 

The techniques currently used in simulations are extremely similar, even 
when the spacetime is evolved using spectral methods ( [^[M] !. multiple patches 
([97]) or a mesh refined cartesian grid (e.g. [98l[99lfT00] ) . In nearly every case a 
finite volume approach combined with a reconstruction-evolution method and 
approximate Riemann solvers is used, computing a conservative update over 
the neutron star and the surrounding artificial atmosphere. These methods 
converge at best at second order and are computationally expensive com- 
pared to the methods required for the spacetime. A notable exception is the 
smoothed-particle hydrodynamics simulations of e.g. [85 11101^102] . but as yet 
these have not used full general relativity (but note the recent results of [103] ). 

These methods have been used to study purely hydrodynamical scenarios. 
The numerical accuracy of the resulting GW signals is not always easy to 
assess, with some work at lower resolutions not claiming accuracy of better 
than 1 orbit ([91]), whilst others at the highest currently available resolutions 
claiming accuracy in key quantities of the order of a percent ( |104| ). Even in 
simulations with careful control of the accuracy in the signal the amplitude 
accuracy is nearly an order of magnitude worse than early black hole binary 
simulations (compare figure 5 of |104j with figure 3 of |105| ). It is also clear 
that numerical accuracy is significantly degraded when shocks appear such as 
in the complex post-merger stage of a binary inspiral. Thus whilst it seems 
that current techniques and computing power are sufficient to make qualitative 
statements about gravitational waves and the general effects of changes in the 
parameter space, detailed quantitative predictions at the level required by 
ET are unlikely to result simply by throwing computing power at current 
techniques, except in certain simple cases. 

Improvements in numerical techniques to significantly improve the accu- 
racy of the results, particularly in the complex post-merger regime, are there- 
fore needed to get the most out of numerical relativity for ET. Extensions of 
standard finite volume techniques to higher order (e.g. |106| ) are one obvious 
avenue to explore. The implementation of more complex but flexible methods 
such as Discontinuous Galerkin finite elements (as implemented in relativity 
by |107pi08) ) is another possibility. Finite difference methods have received 
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some attention (e.g. [lOQpilOj ) and may be the simplest approach. At present 
it seems unhkely that the highly accurate spectral methods will provide a 
solution due to the problems these methods encounter with discontinuities 
(although [111] presents one approach in a particular approximation), but hy- 
brid methods (e.g. |112[|113j l may allow for a combination of the best features 
of spectral or simple high order schemes with HRSC methods. 

Finally it should be noted that the impact of the artificial atmosphere 
used outside the compact GW source on the accuracy of the results remains 
unclear, and methods to avoid the use of an atmosphere (e.g. [114[|115j ) are 
at present not sufficiently advantageous to use in relevant simulations. The 
presence of a free surface (in the perfect ffuid modei) may be an impediment 
to the use of more accurate numerical methods, and remains an area where 
more understanding at both numerical and theoretical level is needed. Purely 
on numerical grounds the use of finite elements is attractive here, as the grid 
can be adapted to the surface of the objects; all that is required is a consistent 
boundary condition to impose. 



3.2 Magnetic fields 

So far the ideal MHD limit has been the focus of a range of studies (for a 
recent review see [IH]), many simply "adding on" the magnetic field to study 
the qualitative differences (see e.g. |1161I117) ). but some studying additional 
instabilities that may arise (e.g. [102illl8j ). With the increase in the parameter 
space due to the strength and topology of the magnetic field the choice of initial 
conditions becomes ever more important, and this is where one fundamental 
difhculty remains. It is expected that the initial object will be approximately 
axi-stationary and the precise topology and mixture of toroidal and poloidal 
components that would make up this initial field has only been studied nu- 
merically under certain assumptions, see [119pi201ll21j . Until a self-consistent 
quasi-stationary solution can be constructed the standard method is to add a 
purely poloidal magnetic field to the interior; it has been argued ( |116| ) that 
this is sufficient for the inspiral phase of a magnetized binary. 

In addition to the conservation constraints, numerical methods should iden- 
tically preserve the V ■ B = constraint. There are many possible approaches, 
with most codes implementing a variant of constrained transport (e.g. [122] ) 
and some using hyperbolic divergence cleaning (e.g. [110| ). Constrained trans- 
port methods should ensure the constraint is satisfied to machine precision 
whilst divergence cleaning propagates any errors away; however divergence 
cleaning is simpler to extend to higher order methods and the presence of 
the neutron star surface and artificial atmosphere may affect the accuracy of 
constrained transport. In all cases the additional constraints to be satisfied 
and fields to be evolved lead to lower numerical accuracy; for example |116j 
shows the numerical error increasing by factors ~ 5 for evolutions of magnetic 
binaries compared to non-magnetic NSs. 
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3.3 Additional effects 

fn the context of computing GWs for ET tlie importance of physics beyond 
the ideal MHD model are the changes to the internal structure, which will 
certainly affect the fine detail and can change the dynamics, leading to a qual- 
itative change in the signal. The most obvious problems arise through thermal 
effects (radiation transport and diffusion post-merger or in core-collapse) or 
where the ideal fluid model breaks down (viscous or multifluid effects such as 
superfluidity, and in the crust). 

3.3.1 Radiation transport 

The transfer of heat and energy through radiation transport changes the local 
structure of the fluid, and is known to be crucial in the dynamics of core- 
collapse supernovae and may play a role in the post-merger dynamics of a 
binary merger. However the simulation of the full Boltzmann equation in 3 -f 
3-1-1 dimensional general relativity would seem to be impractical on timescales 
relevant for ET. An estimate ( |123] ') suggests that even approximate transport 
methods require peta scale computing power. A suggestion for modelling the 
full problem is given by |124| , but on the timescale for ET it seems likely that 
approximation methods such as [125) will be the only practical solution. 

3.3.2 Non-ideal effects 

The consideration of non-ideal effects may become important post-merger 
when the matter surrounding the remnant can be both very hot and yet not 
too dense, leading to a plasma where diffusion and resistivity are important. 
Viscous effects have barely been touched on - for an isolated example see |126j . 
Resistivity effects have been considered and modelled in the magnetosphere, 
but modeUing for neutron stars is only just starting (see |107pi08j ). 

The Newtonian limit gives mixed hyperbolic-parabolic systems of equations 
with their associated causality problems. The relativistic equivalents have stiff 
source terms leading to severe timestep constraints when using standard nu- 
merical methods. Various numerical schemes have been proposed to bypass 
this issue, including the use of IMEX time integrators f [107| V It seems plausi- 
ble that similar numerical techniques will be adopted as those for large mass 
ratio binary black holes. 

3.3.3 Multifluid and elastic effects 

It is believed that most neutron stars have a region containing superfluid 
neutron pairs and that magnetized neutron stars may be superconducting in 
some region. The crust of a mature neutron star is not a fluid but a lattice, 
best modelled as a relativistic elastic solid inter-penetrated by an additional 
fluid. In addition, the possible existence of exotic phases of matter in the inner 
core leads to another region where multiple particles may be flowing freely and 
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independently. All of these effects involve the local interaction of multiple fluids 
or solids (for a review see |127j ). with the appearance of additional dynamics 
and instabilities. 

At present no numerical evolutions of any of these effects have been per- 
formed. Frameworks for evolving multiple fluids have been worked out ( |127| ) 
but have not been tailored for numerical work that would include shocks. It 
seems unlikely that multifluid effects will be visible in numerical simulations 
in the inspiral phase, but the additional propagation speeds and instabilities 
could be visible near or post merger. 



4 Discussion: the nature of numerical relativity in the ET era 

Over the last four decades, numerical relativists have been primarily concerned 
with the technical challenges of solving Einstein's equations numerically: for- 
mulating Einstein's equations to be numerically stable, developing appropriate 
gauge conditions, exploring the numerical techniques necessary to accurately 
evolve relativistic spacetimes, and to handle rclativistic fluids and magnetic 
fields. Many of these issues have been resolved to a sufficient degree that we 
can now extract useful physics from the simulations. In addition, sensitive GW 
detectors arc now in operation and the analysis of their data requires numer- 
ical results. The field of numerical relativity is thus shifting from a focus on 
theoretical and numerical-analysis issues to questions of astrophysics and GW 
data analysis. 

Over the next decade, this shift will continue, and numerical relativists will 
interact more closely with astronomers, astrophysicists, and GW data analysts. 
Closer collaboration between the NR and DA communities has already begun 
through the Numerical INJection Analysis (NINJA) project, where numerical 
black-hole-binary waveforms were injected into simulated Initial LIGO and 
Virgo detector noise, and then the data were analyzed by a battery of search 
and parameter estimation methods, to test their efficacy in dealing with "real" 
GW signals. The successful first NINJA project [128] is now being followed 
with projects to more systematically test search pipelines. Among these are 
the suggestion to use actual detector data, which contain the full array of real 
detector artifacts and can allow the most conclusive tests of search pipelines, 
and to extend the original simulated-dctector-noisc NINJA strategy to matter 
sources (NS-NS, NS-BH and core collapse). 

With the advent of Advanced LIGO/Virgo and the space-based detector 
LISA [129j . the interaction between the NR and DA communities will surely 
become closer, and the possibility of multi-messenger GW astronomy will in- 
clude astronomers and astrophysicists in these efforts, too. Once GW detec- 
tions become routine, and numerical simulations computationally cheaper and 
faster, it may be typical for a candidate observation to motivate follow-up nu- 
merical simulations, which in turn lead to more sophisticated data analysis 
parameter estimation exercises, and comparison with electromagnetic obser- 
vations. The time scale of this back-and-forth interaction may be as short as 
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weeks or even days. For some regions of the (at least black-hole-binary) pa- 
rameter space, it is in fact likely that GW astronomers will have ready access 
to numerical codes that allow them to "dial up" any waveform they require 
for a particular data-analysis exercise, and the focus of numerical-relativity re- 
search will be only on those cases (for example, high mass ratios, binaries with 
complex spin interactions, and matter systems) that still present problems. 

For matter simulations considerable work has been done to explore the 
parameter space by including as much of the relevant physics as is currently 
practical. However, in order to take full advantage of the additional accuracy 
given by ET these simulations will have to improve substantially in certain 
areas. Firstly the numerical methods currently employed will probably not 
deliver the improvements in accuracy required based on the (conservative) 
increase in computing power projected here. Secondly the community will have 
to collaborate closely with the data analysis community in order to determine 
which aspects of the physical model will actually lead to interesting observable 
effects, and hence where in the parameter space the numerical simulations 
should be focused. As the discussion in section [3] makes clear it would be 
possible to spend all the time and increased computing power in applying and 
simulating better physical models and exploring the parameter space. It seems 
likely that to take the greatest advantage of ET a narrower focus is needed to 
produce the most accurate simulations including the most relevant physics for 
the best candidate sources. 
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